Direct Sum

A direct sum collects several blocks into one state. This is useful when the unknowns have different meanings or different tensor spaces, such as a symmetric stress tensor and a scalar internal variable.

In Tensorial, direct-sum values are built with pack. The infix operator ⊕ is equivalent to pack.

Pack and unpack blocks

julia> A = @Mat [1.0 2.0; 3.0 4.0]2×2 Tensor{Tuple{2, 2}, Float64, 2, 4}:
 1.0  2.0
 3.0  4.0
julia> x = pack(A, 3.0)2-element DirectSumVector with storage Float64:
 Space(2, 2)
 Space()
julia> unpack(x)([1.0 2.0; 3.0 4.0], 3.0)
julia> unpack(x, 1)2×2 Tensor{Tuple{2, 2}, Float64, 2, 4}:
 1.0  2.0
 3.0  4.0
julia> unpack(x, 2)3.0

The infix form is convenient for short expressions:

julia> y = A ⊕ 3.02-element DirectSumVector with storage Float64:
 Space(2, 2)
 Space()
julia> x == ytrue

For symmetric tensor blocks, the flat coordinate representation uses Mandel scaling. This is visible with flatview:

julia> As = symmetric(A)2×2 SymmetricSecondOrderTensor{2, Float64, 3}:
 1.0  2.5
 2.5  4.0
julia> z = pack(As, 3.0)2-element DirectSumVector with storage Float64:
 Space(Symmetry(2, 2),)
 Space()
julia> flatview(z)4-element StaticArraysCore.SVector{4, Float64} with indices SOneTo(4):
 1.0
 3.5355339059327378
 4.0
 3.0

Differentiating packed states

Several coupled variables can be treated as one state while derivatives keep the block layout. Here the state contains a symmetric tensor block A and a scalar block s:

julia> x = pack(symmetric(@Mat [1.0 2.0; 2.0 4.0]), 2.0)2-element DirectSumVector with storage Float64:
 Space(Symmetry(2, 2),)
 Space()
julia> function f(z)
           A, s = unpack(z)
           dot(A, A) + s * tr(A) + s^2
       endf (generic function with 1 method)
julia> f(x)39.0

The gradient has the same block structure as x:

julia> g = gradient(f, x)2-element DirectSumVector with storage Float64:
 Space(Symmetry(2, 2),)
 Space()
julia> unpack(g)([4.0 4.0; 4.0 10.0], 9.0)

The Hessian is a direct-sum matrix. The tensor–tensor block is fourth-order, so here we check its type and print the smaller coupling blocks:

julia> H = hessian(f, x)2×2 DirectSumMatrix with storage Float64:
 Space(Symmetry(2, 2), Symmetry(2, 2))  Space(Symmetry(2, 2),)
 Space(Symmetry(2, 2),)                 Space()
julia> unpack(H, 1, 1) isa SymmetricFourthOrderTensor{2}true
julia> unpack(H, 1, 2)2×2 SymmetricSecondOrderTensor{2, Float64, 3}:
 1.0  0.0
 0.0  1.0
julia> unpack(H, 2, 1)2×2 SymmetricSecondOrderTensor{2, Float64, 3}:
 1.0  0.0
 0.0  1.0
julia> unpack(H, 2, 2)2.0

Residuals and Jacobian blocks

Packed states also work for residual maps, where the output is packed as well:

julia> C = symmetric(@Mat [2.0 1.0; 1.0 3.0])2×2 SymmetricSecondOrderTensor{2, Float64, 3}:
 2.0  1.0
 1.0  3.0
julia> function F(z)
           A, s = unpack(z)
           B = A + s * C
           t = tr(A) + s^2
           pack(B, t)
       endF (generic function with 1 method)
julia> r = F(x)2-element DirectSumVector with storage Float64:
 Space(Symmetry(2, 2),)
 Space()
julia> unpack(r)([5.0 4.0; 4.0 10.0], 9.0)

With

\[\bm{z} = \begin{bmatrix} \bm{A} \\ s \end{bmatrix}, \qquad \bm{F}(\bm{z}) = \begin{bmatrix} \bm{B}(\bm{A},s) \\ t(\bm{A},s) \end{bmatrix},\]

the Jacobian has the corresponding block structure:

\[\bm{J} = \frac{\partial \bm{F}}{\partial \bm{z}} = \begin{bmatrix} \dfrac{\partial \bm{B}}{\partial \bm{A}} & \dfrac{\partial \bm{B}}{\partial s} \\ \dfrac{\partial t}{\partial \bm{A}} & \dfrac{\partial t}{\partial s} \end{bmatrix}.\]

The calls below inspect the same four Jacobian blocks. Again, the first block is a fourth-order tensor, while the remaining blocks are small enough to print:

julia> J = gradient(F, x)2×2 DirectSumMatrix with storage Float64:
 Space(Symmetry(2, 2), Symmetry(2, 2))  Space(Symmetry(2, 2),)
 Space(Symmetry(2, 2),)                 Space()
julia> unpack(J, 1, 1) isa SymmetricFourthOrderTensor{2}true
julia> unpack(J, 1, 2)2×2 SymmetricSecondOrderTensor{2, Float64, 3}:
 2.0  1.0
 1.0  3.0
julia> unpack(J, 2, 1)2×2 SymmetricSecondOrderTensor{2, Float64, 3}:
 1.0  0.0
 0.0  1.0
julia> unpack(J, 2, 2)4.0

This is the structure needed by Newton updates and other coupled local solves. If J is a DirectSumMatrix and r is a DirectSumVector, the correction can be written directly:

julia> δx = -J \ r2-element DirectSumVector with storage Float64:
 Space(Symmetry(2, 2),)
 Space()
julia> unpack(δx)([7.0 2.0000000000000004; 2.0000000000000004 8.0], -6.0)

No manual flattening is needed.

julia> flatview(J)4×4 StaticArraysCore.SMatrix{4, 4, Float64, 16} with indices SOneTo(4)×SOneTo(4):
 1.0  0.0  0.0  2.0
 0.0  1.0  0.0  1.41421
 0.0  0.0  1.0  3.0
 1.0  0.0  1.0  4.0

Example: return mapping residual

As a larger example, consider the active plastic branch of a small-strain isotropic J2 return-mapping update. The unknown state contains

• the updated symmetric stress σ, and
• the plastic multiplier increment Δγ.

We solve the local residual

\[\bm{R}(\bm{\sigma}, \Delta\gamma) = \begin{Bmatrix} \bm{\sigma} - \bm{\sigma}^{\mathrm{tr}} + \Delta\gamma\,\bm{\mathbb{C}}^{\mathrm{e}} : \bm{n} \\ q(\bm{\sigma}) - (\sigma_{y0} + H\,\Delta\gamma) \end{Bmatrix} = \bm{0}.\]

Here σᵗʳ is the trial stress, ℂᵉ is the elastic stiffness, q is the von Mises stress, σy0 is the initial yield stress, and H is the isotropic hardening modulus. The flow direction n is the stress derivative of the yield function.

The chosen trial stress is in the plastic branch. If q(σᵗʳ) ≤ σy0, the update is elastic: σ = σᵗʳ and Δγ = 0.

σᵗʳ = SymmetricSecondOrderTensor{3}((2.0, 0.4, 0.2, 1.2, 0.1, 0.9))
K = 10.0  # bulk modulus
G = 5.0   # shear modulus
σy0 = 0.6 # initial yield stress
H = 2.0   # isotropic hardening modulus

Ivol = vol(SymmetricFourthOrderTensor{3}) # volumetric projector
Idev = dev(SymmetricFourthOrderTensor{3}) # deviatoric projector
ℂᵉ = 3K * Ivol + 2G * Idev

q(σ) = sqrt(3/2) * norm(dev(σ)) # von Mises stress
yield_function(σ, Δγ) = q(σ) - (σy0 + H * Δγ)

function R(x)
    σ, Δγ = unpack(x)
    # flow direction and yield-function value
    n, f = gradient(σ -> yield_function(σ, Δγ), σ, :all)
    R_σ = σ - σᵗʳ + Δγ * (ℂᵉ ⊡₂ n)
    pack(R_σ, f)
end

x = pack(σᵗʳ, 0.0)
unpack(R(x))

([0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0], 0.6649110640673516)

The Jacobian is obtained directly from the packed residual:

J = gradient(R, x)
n = gradient(σ -> yield_function(σ, 0.0), σᵗʳ)

(
    unpack(J, 1, 1) ≈ one(SymmetricFourthOrderTensor{3}),
    unpack(J, 1, 2) ≈ ℂᵉ ⊡₂ n,
    unpack(J, 2, 1) ≈ n,
    unpack(J, 2, 2) ≈ -H,
)

(true, true, true, true)

A Newton correction can then be written in direct-sum form:

δx = -J \ R(x)
xnew = x + δx
σ, Δγ = unpack(xnew)
norm(R(σ ⊕ Δγ))

4.263892432392673e-16

For this radial-return example, one Newton update gives the return-mapping solution, so the residual is zero for the updated state. More general return-mapping problems may need several Newton iterations.

For direct-sum docstrings, see Direct sum API.